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ABSTRACT 


Dynamics  of  the  circulation  in  the  Sea  of  Marmara  are 
investigated  with  a time  dependent,  three  dimensional 
numerical  model.  The  empirically  inferred  hydrologic  regimes 
of  the  sea  and  connective  straits  are  discussed. 

A baroclinic  model  based  on  the  primitive  equations  is 
solved  by  direct  integration  of  an  initial  value  problem. 

The  circulation  in  the  Sea  is  driven  by  surface  forces  that 
simulate  wind  stress  and  horizontal  pressure  gradient  forces 
related  to  internal  stratification. 

The  predicted  density  field,  in  sigma— t units, is  compared 
with  data.  Detailed  three  dimensional  horizontal  velocity 
patterns  and  vertical  velocity  patterns  in  horizontal  planes 
are  given. 

Bottom  friction,  irregular  bottom  topography,  non— linear 
terms  in  the  momentum  equation  and  vertical  mean  part  of  the 
horizontal  velocity  have  been  omitted.  For  simplicity 
density  is  predicted  in  place  of  temperature  and  salinity. 
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I . INTRODUCTION 


There  is  an  intense  growing  interest  in  recent  years  in 
the  development  of  simulation  of  hydrological  conditions  in 
coastal  waters  and  smaller  adjacent  seas.  Advancement  in 
modern  technology  requires  much  greater  accuracy  than  can  be 
obtained  by  the  use  of  classical  instrumentation  and  measure- 
ment techniques  only.  Due  to  the  discontinuous  nature  of 
such  observations  these  measurements  can  provide  the  basis 
for  analysis  of  the  current  system,  temperature,  salinity 
and  other  variables  on  a climatological  basis.  But  climato- 
logical conditions  are  not  accurate  enough  in  coastal  waters 
for  uses  such  as  construction,  navigation,  tracing  of  pollu- 
tants, mining  effluent,  oil  spills,  search  and  rescue  opera- 
tions, agriculture,  etc.  It  is  also  true  that  sampling  the 
ocean  adequately  whether  for  exploratory  or  research  purpose 
is  a major  problem.  Great  short-term  variability  of  condi- 
tions in  coastal  waters  makes  it  even  more  difficult  to  make 
adequate  representative  measurements. 

As  a result,  in  addition  to  the  classical  methods  of 
analyses  of  current  systems,  determination  of  salinity  and 
temperature  distributions  etc.,  new  simulation  techniques 
which  were  applied  long  before  synoptic  numerical  weather 
prediction  received  much  more  attention  (T.  Levastu,  S. 

Larson  and  K.  Rabe,  1976) . 
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The  objective  of  this  study  is  to  develop  a regional 
circulation  model  of  the  Sea  of  Marmara.  The  circulation 
and  density  distribution,  in  sigma— t units,  in  that  sea  are 
simulated  by  a nine— level  numerical  model  with  an  hori- 
zontal flat  bottom  at  1000  m depth.  The  motion  is  driven 
by  prescribed  wind  stress  and  horizontal  pressure  gradient 
forces  due  to  internal  stratification. 


ll 


II.  PHYSICAL  GEOGRAPHY  AND  HYDROLOGY  OF 


THE  SEA  OF  MARMARA 

The  Sea  of  Marmara  is  a small  sea  located  between  Asia 
Minor  and  Europe  and  occupies  the  main  part  of  the  Turkish 
straits  system  (210  km) . The  Strait  of  Canakkale  (Darda- 
nelles) (60  km)  connects  the  Sea  of  Marmara  to  the  Mediter- 
ranean and  through  the  Strait  of  Istanbul  (Bosporus)  (30  km) 
it  is  connected  to  the  Black  Sea  at  the  North.  (Figure  1) 

The  structure  of  the  basin  which  is  small  in  area  (about 
2 

11000  km  ) is  associated  with  the  Anatolicin  fault  which 
lies  along  its  bottom  to  the  North.  The  trough  passing  along 
the  Northern  deep  slope  consists  of  three  depressions  (max. 
1380  m) . These  three  basins  are  separated  by  low  connecting 
sills.  A shallow  trough  lies  at  the  foot  of  the  continental 
slope  which  is  well  pronounced  only  in  the  easternmost  part 
of  the  sea.  At  the  South  it  becomes  shallower  having  depths 
of  the  order  60—80  m.  Minimum  depths  at  the  Strait  of 
Canakkale  and  Strait  of  Istanbul  are  57  and  37  m;  they  have 
mean  widths  4.5  km  and  0.7  km  and  average  lengths  60  and  30 
km,  respectively. 

The  climate  of  the  Sea  of  Marmara  is  influenced  by  the 
regional  atmospheric  circulation  system  of  the  Eastern  Medi- 
terranean, the  Black  Sea  and  Asia  Minor  Peninsula.  During 
the  summer,  the  Eastern  Mediterranean  basin  is  under  the 


Figure  1.  Location  map  of  the  Sea  of  Marmara 


influence  of  a stable  wind  system  which  is  called  "Etesian" 
wind.  Etesians  are  monsoon— like  and  associated  with  a quasi— 
stationary  cyclonic  circulation  above  the  Asia  Minor  Penin- 
sula during  summer.  The  active  period  of  the  Etesian  winds 
falls  almost  entirely  within  the  summer  season  (approximate- 
ly from  May  to  October) . There  are  two  maximums,  which  are 
pronounced  during  July  and  August.  These  Etesian  winds, 
which  are  dry  winds  and  generally  northerly  or  northeasterly 
from  the  Black  Sea,  have  diurnal  variability. 

Southerly  winds  bring  warm  and  humid  air  to  the  area 
during  winter.  When  these  southerly  winds  are  sufficiently 
strong,  they  produce  storm  waves  and  surges  along  the  coast 
of  Marmara  and  at  the  entrance  to  the  Strait  of  Istanbul 
(Gunnerson  and  Ozturgut,  1974) . 

The  hydrological  regime  of  the  Sea  of  Marmara  is  largely 
determined  by  the  water  exchange  between  the  Mediterranean 
and  Black  Seas  through  the  connective  straits.  Strait  of 
->  Canakkale,  and  Strait  of  Istanbul.  Problems  of  water  ex- 
change, vertical  mixing  and  hydrological  characteristics  of 
the  currents  in  these  straits  and  in  the  Sea  of  Marmara  have 
been  the  subject  of  many  extensive  studies  conducted  by  the 
Turkish  Navy  Hydrographic  Office  and  by  the  University  of 
Istanbul.  In  addition  to  these  fundamental  and  frequent 
studies  there  were  some  investigations  conducted  by  German, 
Russian  and  English  investigators  at  the  end  of  the  19th 
and  beginning  of  the  20th  century.  Recent  studies  and  obser- 
vations, if  connected  to  old  data  and  observations,  can  give 
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a better  understanding  of  the  complete  regional  oceanographic 
condition. 

These  studies  are  mainly  concentrated  on  the  Strait  of 
Istanbul  and  the  vicinity  of  this  strait.  This  is  due  to  the 
importance  of  this  strait  as  a link  in  the  connection  between 
the  Black  Sea  and  the  Mediterranean  Sea  through  the  Sea  of 
Marmara.  This  makes  it  an  important  factor  controlling  Black 
Sea  chemistry,  biology,  sedimentation,  etc. 

-*>  Based  on  the  first  detailed  survey  of  Makarov  and  Merz 

1917—18,  in  Moller  (1928)  , and  subsequent  studies  of  Ullyott 
->  and  Ilgaz  (1946),  Pektas  (1956),  A.  K.  Bogdanova  (1961,  1965), 
C.  G.  Gunnerson  and  E.  Ozturgut  (1974)  and  others,  the  main 


properties  of  this  straits'  system  and  its  vicinities  are 
fairly  well  understood. 

One  of  the  main  characteristics  of  these  straits  is  a 
two— layer  structure  of  the  current,  which  is  determined  by 
the  differences  between  the  relatively  fresh  water  of  the 
Black  Sea  and  saline  water  from  the  Mediterranean  which  occu- 
pies the  depths  of  the  Sea  of  Marmara  and  forms  the  bottom 
water  of  two  connective  straits  (Figure  2) . Other  external 
factors  such  as  wind,  water  level  drop  along  the  straits, 
water  balance,  topographic  influences  and  the  depth  and  width 
of  the  straits  also  must  be  considered.  The  scale  of  the 
motion  and  complexity  of  these  factors  make  it  difficult  to 
determine  the  current  system  mathematically. 

Another  main  characteristic  of  these  straits  is  the  two- 
layer  stratification  of  water  which  is  associated  with  the 
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two-layer  current  system.  This  stratification  is  also  a 
consequence  of  the  water  mass  differences  at  the  two  ends 
of  the  straits. 

Detailed  observations  of  temperature  and  salinity  dis- 
tributions indicate  certain  features  of  the  diffusion  of 
the  Mediterranean  water  into  the  Black  Sea.  This  water 
flows  North  through  the  Strait  of  Istanbul  from  the  Sea  of 
Marmara.  It  is  shown  by  Bogdanova  (1961,  1965)  that  this 
underflow  of  Mediterranean  water  into  the  Black  Sea  is 
present  throughout  the  year  and  descends  to  significant 
depths  in  the  Black  Sea. 

This  water  exchange  through  the  straits  system  fluctu- 
ates according  to  seasons,  years  and  periods.  This  is  very 
important  for  the  hydrologic  regimes  of  the  Sea  of  Marmara 
and  the  Black  Sea.  Water  level  differences,  given  with  the 
following  expression,  are  mainly  due  to  excess  precipita- 
tion over  evaporation  and  river  runoff  in  the  Black  Sea. 


ZB  ZM 

sea  level  differential  between 
Black  Sea  and  Mediterranean 

water  level  of  Black  Sea 

water  level  of  Mediterranean 


(1.1) 


Sea  level  differential  (Ah)  is  always  positive  and  fluctu- 
ates with  time.  It  can  be  expressed  in  two  components 


Ah  = AK  + (Ah) ' 


(1.2) 
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where 


ATT  : average  sea  level  differential 

(Ah) 1 : fluctuating  part  of  sea  level 
differential  [(Ah)'  « f(t)] 

According  to  Bogdanova  (1965)  Ah  has  a value  42  cm 
based  on  the  long  term  observations.  The  fluctuating  part 
of  the  sea  level  differential  and  total  level  differential 
are  given  in  Table  I. 


TABLE  I. 

Monthly  average  fluctuation  sea  level  differential 
and  sea  level  differential.! 


Month 

I 

II 

III 

IV 

V 

VI 

VII 

VIII 

IX 

X 

XI 

XII 

(Ah) ' (cm) 

1 

0 

3 

7 

11 

15 

6 

0 

-5 

-7 

-7 

-4 

Ah  (cm) 

43 

42 

45 

49 

53 

57 

48 

42 

37 

35 

35 

38 

(^Based  on  the  data  given  by  Bogdanova,  1965) 


The  maximum  change  is 


8 1 


(4hlmax  = Tt  (ih')«,ax  “ “ 9 cm/inonth' 


observed  during  approximately  June.  This  falls  approximately 
in  a period  when  continental  runoff  begins  to  decrease.  The 
fluctuating  sea  level  differential  increases  from  March  to 
July.  As  a consequence  during  this  period  the  upper  current 
becomes  stronger  and  the  lower  current  becomes  weaker.  From 
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August  to  November  sea  level  fluctuation  decreases  and  be- 
gins to  increase  again  from  August  to  February.  The  upper 
current  of  the  Strait  of  Istanbul  becomes  weaker  and  the 
lower  current  becomes  stronger. 

According  to  Bogdanova's  (1965)  observations  the  maxi- 
mum Ah  value  is  observed  during  June  and  has  a value  57 
cm;  the  minimum  Ah  value  35  cm  is  observed  during  October 
and  November. 

The  average  sea  level  differential  between  the  Black 
Sea  and  the  Sea  of  Marmara  is  about  35  cm  which  is  six 
times  Moller's  (1928)  estimate  (Gunnerson  and  Ozturgut, 

1974)  . 

There  is  a characteristic  stratification  with  a sharp 
density  transition  layer  over  the  whole  area.  The  sharpness 
is  more  pronounced  at  the  North  at  the  Strait  of  Istanbul. 

A cross  section  of  sigma— t values  for  the  Strait  of  Istanbul 
is  shown  in  Figure  3.  Fresh  warm  water  of  the  Black  Sea 
lies  in  a surface  layer  15—20  m thick  at  Istanbul.  This 
water  which  is  mainly  river  runoff  and  excess  precipitation 
overrides  the  Mediterranean  water  which  is  cold  and  dense. 

The  wedge  form  of  the  upper  water  shows  clearly  in  this 
strait,  but  this  is  not  true  for  the  Strait  of  Canakkale. 

On  the  other  hand,  the  wedge  form  of  the  lower  water  is  com- 
paratively stronger  and  more  pronounced  in  the  Strait  of 
Canakkale  due  to  the  downward  sloping  of  the  sea  bed  toward 
the  North  in  the  Strait  of  Istanbul.  Average  sigma— t values 
for  surface  and  bottom  water  are  13.5  and  27.5.  (Defant,  1961) 
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According  to  Gunnerson  and  Ozturgut  (1974) , the  average 
salinity  of  Black  Sea  waters  at  a station  9 km  from  the 
Northern  entrance  of  the  strait  of  Istanbul  was  approximate- 
ly 17.5  °/oo  for  an  observation  period  July  1966  to  December 
1967.  Lower  values  were  observed  (16—17  °/oo)  from  July  17 
to  August  23,  1967.  This  lower  surface  salinity  reflects 
maximum  discharge  of  rivers  into  the  Black  Sea.  Higher 
average  surface  salinities  (17.5  — 19.0  °/oo)  were  observed 
for  the  rest  of  the  year  and  extreme  values  (as  high  as 
25  °/oo)  during  winter  months.  These  were  due  to  the  favor- 
able conditions  for  Mediterranean  water  to  flow  northward. 
These  conditions  include  small  differences  between  the  levels 
of  the  Black  Sea  and  the  Mediterranean  and  southerly  or 
southwesterly  winds  over  the  Sea  of  Marmara  and  Strait  of 
Istanbul.  At  the  midpoint  the  average  salinity  in  the  lower 
layer  is  approximately  38.5  °/oo  and  in  the  surface  layer 
17.5  °/oo.  A salinity  cross  section  is  given  in  Figure  4. 

Warm  and  fresh  Black  Sea  water  can  be  identified  easily 
from  the  temperature  cross-section  which  is  given  in  Figure 
5.  Below  this  warm  fresh  water  of  the  Black  Sea  cold  inter- 
mediate water  of  Sea  of  Marmara  can  be  observed.  Mediter- 
ranean water  occupies  the  deep  bottom  layer  of  the  Strait  of 
Istanbul. 

The  dynamics  of  the  processes  occurring  in  these  two 
straits  is  determined  by  the  wind  stress  over  the  straits 
and  the  drop  of  the  sea  levels  and  densities  at  the  end  of 
the  straits.  The  main  driving  forces  are  gradient  of  water 
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level,  thermohaline  forcing,  wind  force  exerted  on  the  sur- 
face and  coriolis  force. 

According  to  observations  the  boundary  between  the  water 
masses  generally  does  not  coincide  with  the  boundary  between 
currents.  The  boundary  between  water  masses  has  a greater 
slope  than  the  boundary  between  currents  (Defant,  1961)  . 

This  is  probably  due  to  diffusion  processes  occurring  at  the 
current  boundary  and  to  the  nature  of  the  current  (Figure  2) . 

The  average  Black  Sea  inflow  and  outflow  through  the 
Straits  of  Istanbul  according  to  the  analysis  of  Moller  (1928) 
are  6,100  and  12,600  m^/sec,  respectively  (Gunnerson  and 
Ozturgut,  1974) . Velocity  is  greatest  at  the  sea  surface 
and  decreases  rapidly  with  depth  and  laterally;  it  increases 
from  North  to  South.  Under  normal  conditions  it  is  40—50 
cm/sec  at  the  entrance  of  the  strait  and  150  cm/sec  at  the 
South  end.  The  lower  current  is  strongest  in  the  central 
part  of  the  lower  water  which  is  at  about  16  m from  the 
bottom  in  the  Strait  of  Istanbul  and  45  m in  the  Strait  of 
Canakkale.  These  velocities  are  100—150  cm/sec  and  25  to 
10  cm/sec  respectively  (Defant,  1961) . 

In  general  in  the  Sea  of  Marmara  there  is  Black  Sea 
water  at  the  surface  and  Mediterranean  water  at  the  bottom. 
This  is  the  normal  state.  Except  for  times  when  very  strong 
southerly  winds  bring  bottom  water  to  the  surface  and  cause 
wind  mixing  these  two  main  water  types  are  both  present  with 
a very  thin  intermediate  water  layer.  The  Black  Sea  water 
extends  to  a depth  of  15—20  m in  the  vicinity  of  the  Strait 
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of  Istanbul  and  becomes  shallower  farther  south.  There  is 
intermediate  water  between  these  two  water  types  which  ex- 
tends at  most  to  about  50—60  m.  Below  this  the  basin  is 
filled  with  Mediterranean  water  as  shown  in  Figures  6,7  and 
8.  This  water  has  a temperature  14.9°C,  salinity  38.55  °/oo 
and  sigma— t 28.75,  approximately  in  June. 

Surface  water  has  a large  annual  variability.  On  the 
other  hand  bottom  water  does  not.  Due  to  the  described 
oceanographic  conditions  bottom  water  is  almost  isolated  from 
surface  processes  and  has  a large  time— scale  of  variability 
compared  to  surface  water.  Average  inflow  and  outflow  rates 
at  the  straits  give  approximately  30  years  for  "turnover" 
time  or  "flushing"  time  for  bottom  water  and  200  days  for 
surface  water. 

Based  on  the  same  rates  a typical  value  for  the  vertical 
mean  current  |u|  is  10  cm  sec  . A characteristic  magni- 
tude for  horizontal  velocity  |u|  using  the  thermal  wind 
relation  is  4 cm  sec  ^ . Due  to  the  large  isolated  reservoir 
of  Mediterranean  water  of  almost  uniform  characteristics  in 
the  lower  part  of  the  basin,  the  thermohaline  circulation 
is  confined  to  the  upper  50—60  m. 
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Figure  6.  Sigma— t cross  section  of  Sea  of  Marmara  for  upper 
100  meters.  Section  is  taken  between  Strait  of  Canakkale  an' 
Strait  of  Istanbul  (approximately  B— C point  in  figure  2) . 


Figure  7.  Temperature  cross  section  of  Sea  of  Marmara  for  upper 
100  meters.  (Same  section  as  Figure  6.) 


III.  SELECTION  OF  THE  NUMERICAL  MODEL 


Regional  circulation  models  may  be  separated  into  two 
categories.  One  type  of  regional  circulation  model  is  a 
diagnostic  density  model.  This  approach  is  based  on  an  ob- 
served density  field  and  wind  stress  distribution.  These 
quantities,  however,  are  not  observed  continuously  and  in 
an  evenly  distributed  way.  But  based  on  average  values  of 
these  quantities  a velocity  field  is  calculated  diagnosti- 
cally. A diagnostic  model,  which  seems  at  first  easier, 
has  drawbacks  mainly  of  two  types.  First,  it  requires  tak- 
ing into  account  actual  configuration  of  the  basin  geometry 
and,  second,  it  requires  sufficient  representative  observa- 
tions. These  two  factors  make  this  approach  complicated 
and  expensive. 

A second  type  of  approach  to  study  regional  circulation 
models  is  to  use  a predictive  model  with  an  idealized  topo- 
graphy and  basin  area.  But  this  becomes  more  complex  by  the 
inclusion  of  more  dynamical  processes. 

At  first  glance  it  seems  that  diagnostic  models  are  a 
better  way  of  simulation  of  dynamical  processes.  But  these 
observations  on  which  such  models  depend  are  generally 
gathered  by  land  based  stations  (e.g.,  wind  data)  and  extra- 
polated to  the  area  of  concern.  Furthermore,  using  average 
values  brings  up  the  question  of  how  representative  the 
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result  is.  Further,  it  is  difficult  to  get  a continuous 
view  of  the  dynamical  process  and  long-term  quantative 
analyses . 

On  the  other  hand,  although  prognostic  models  are  complex 
due  to  the  dynamical  processes  involved,  they  are  less  expen- 
sive and  can  give  a continuous  view  of  variables,  which  is 
generally  based  on  an  initial  density  field.  At  least  as  a 
first  approximation  it  seems  reasonable  and  economical  to 
study  a regional  circulation  simulation.  Present  data  are 
still  used  to  estimate  the  initial  density  field  and  wind 
stress  distribution  in  this  kind  of  model. 

Due  to  the  uneven  distribution  of  density  and  wind  stress 
observations  it  was  mandatory  to  use  the  prognostic  approach 
at  the  beginning.  The  procedure  adopted  is  to  use  present 
data  as  efficiently  as  possible  to  estimate  a realistic  ini- 
tial density  field  and  a constant  wind  stress  distribution 
throughout  the  integration  period. 


IV.  DESCRIPTION  OF  THE  DYNAMIC  MODEL 

A.  BASIC  EQUATIONS  OF  THE  MODEL 

Basic  assumptions  which  lead  to  important  simplifications 
are:  (1)  the  hydrostatic  assumption,  (2)  the  Boussinesq  assump- 
tion with  respect  to  density  variations  (density  variation  is 
neglected  except  where  it  appears  as  a coefficient  of  the 
gravitational  constant) , and  (3)  an  implicit  treatment  of  the 
mixing  process  by  eddy  diffusion.  In  addition,  the  non— linear 
terms  are  neglected  in  the  momentum  equation. 

With  these,  the  sea  is  assumed  incompressible,  and  density 
is  replaced  by  a constant  value  everywhere  except  in  the  hydro- 
static equation. 

The  governing  equations  based  on  these  assumptions  are: 
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where  the  symbols  have  the  following  meanings: 
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time 

horizontal  coordinate  measures  positive  eastward 

horizontal  coordinate  measures  positive  northward 

vertical  coordinate  measures  positive  upward  with 
an  origin  at  mean  water  level 

eastward  velocity  component 

northward  velocity  component 

vertical  velocity  component 

density 

reference  constant  density 

departure,  from  reference  pressure  which  is 

pressure  at  reference  density  (=  — pQgz) 

acceleration  of  gravity 

coriolis  parameter  (=  2 ft  sincf)) 

angular  speed  of  earth  rotation 

latitude 

coefficient  of  horizontal  diffusion  for  momentum 
(constant) 

coefficient  of  vertical  diffusion  for  momentum 
(constant) 

coefficient  of  horizontal  diffusion  for  density 
(constant) 

coefficient  of  vertical  diffusion  for  density 
(constant) 

material  derivative  [=  + u + v + w ] 

introduces  the  effect  of  convective  process  when 
the  stratification  is  unstable 
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The  sea  is  driven  by  wind  and  thermohaline  forcing,  these 
effects  being  introduced  with  surface  boundary  conditions  in 
the  momentum  and  density  equations,  respectively.  The  bound- 
ary conditions  specified  at  the  sea  surface  are 


V 7z 


v 3z 


z = 0 


(4.6) 


Heating  or  cooling  at  the  sea  surface  is  neglected  over 
the  period  of  integration.  A linear  north— south  density  gra- 


dient is  specified  in  the  uppermost  level  in  the  model  and 
held  constant  during  calculations.  The  x component  of  sur- 
face stress  is  taken  to  be  zero  and  calculations  are  carried 
out  with  a uniform  meridional  wind  stress. 

The  boundary  conditions  specified  at  the  flat  horizontal 
bottom  are 


9p  _ 


z = -H 


(4.7) 


u=v=w=  0 





The  boundary  conditions  at  the  side  walls  of  the  basin 


are 


^ |£  = u = v = 0 


ah|£'u  = v = 0 


y = o,  Ly 


x = 0 , Lx 


(4.8) 


where  Ly  : width  of  the  basin  ( 50  Jem) 

Lx  : length  of  the  basin (200  Jem) 


Boundary  conditions  at  the  straits  are  based  on  specified 
density  and  velocity  fields  which  were  taken  from  available 
data  and  annual  average  inflow  and  outflow  rates  respectively. 


B.  DESCRIPTION  OF  THE  SOLUTION  TECHNIQUES 

Since  the  model  sea  has  a horizontal  flat  bottom  and 
"rigid  lid"  approximation,  vertical  velocity  vanishes  at  the 
surface  and  bottom.  Following  the  method  used  by  Bryan  and 
Cox  (1967— 1968a— b) , Haney  (1974)  and  others,  the  ocean 
surface  was  considered  a balanced  surface  at  which  w = 0. 

In  this  approach  the  height  of  the  sea  surface  is  not  obtained 
from  the  continuity  equation;  therefore,  external  gravity 
waves  are  eliminated.  This  also  removes  the  vertical  mean 
divergence  from  all  scales  of  motion.  Elimination  of  external 
gravity  waves  allows  the  use  of  a longer  time  step.  In  the 
case  of  long-term  integration  this  is  an  important  gain  com- 
pared to  the  loss  of  the  prediction  for  height  of  the  sea 
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surface.  Filtering  external  gravity  waves  has  a negligible 
effect  on  the  density— driven  part  of  the  current. 

Due  to  the  conditions  on  vertical  velocity,  integration 
of  the  continuity  equation  from  bottom  to  surface  requires 
the  following  condition: 
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where  the  meaning  of  the  overbar  is 
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( ) vertical  mean  [= 


-iS 


( ) dz] 
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and 


( )' indicates  departure  from  vertical  mean 
of  any  quantity. 


Also  for  some  basic  integral  constraints  this  is  a required 
relation  to  be  consistent  with  boundary  conditions. 

The  "rigid  lid"  approximation,  and  (4.9),  are  accomplished 
by  vertically  integrating  equations  (4.1),  (4.?)  and  subtract- 

ing the  result  from  equations  (4.1)  and  (4.2)  respectively. 
Using  present  boundary  conditions,  the  new  sets  of  equations 
are 
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where 


u'  eastward  shear  velocity  component 

v'  northward  shear  velocity  component 

P • departure  of  P from  its  vertical  average 

[=  f (P-P0)gdn  - | f°  ( f°  (p-p0)gdn)dz]  (4.11) 

z —Hz 

The  relation  between  shear  current  and  total  current  for 
horizontal  components  is 

u' (x,y , z , t)  = u (x,y , z , t)  - u(x,y,t) 

(4.12) 

v'(x,y,z,t)  = v(x,y,z,t)  - v(x,y,t) 


since 

| u | > > [ u [ and  ( v [ > > [ v | 
equation  (4.12)  takes  the  form 
u' (x,y,z,t)  = u(x,y,z,t) 

v'(x,y,z,t)  = v (x,y , z , t) 


(4.13) 


(4.14) 


Thus  total  u,  v component  of  horizontal  velocity  can  be 
represented  by  vertical  shear  currents  u',  v'.  In  the 
following  part  primes  of  u'  and  v'  are  dropped. 
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The  equations  are  non— dimensionalized  by  making  the  fol- 
lowing substitutions  for  new  non-dimensional  variables 


(x,y)  = (x,y)/L 
(z)  = (z)/h 

(u,v)  = (u,v)/V. 


(w)  = (wj/Vjh 


where 


/p  p.)= 

« fru 
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(4.15) 


(t) 

L 

V1 

_ (t) 


scale  velocity  associated  with  density  driven 
current 

scale  depth,  taken  equal  to  a characteristic 
thermocline  depth 

length  scale 

X V 

wind  stress  [=  MAX(x  (x,y) , TJ(x,y)] 

3 

sigma— t [=  (p— 1)  x 10  ] 

surface  sigma— t at  northern  part  of  the  basin 
surface  sigma— t at  southern  part  of  the  basin 


■ * . *1 


A scale  velocity  associated  with  the  density  driven  cur- 
rent is  defined  by  using  geostrophic  and  hydrostatic  scaling 
which  were  introduced  by  Bryan  and  Cox  (1968a) . When  the  geo— 
strophic  relation  is  differentiated  with  respect  to  z and  the 
hydrostatic  equation  is  used,  the  result  is 


2.  3P  - f 9u 
PQ  9y  3z 


(4.16) 


A scale  velocity  connected  with  the  density  distribution 
may  be  defined  by 
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where 


Ap 


north— south  surface  density  differences 
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A scale  velocity  for  wind  driven  current  is  defined  by 


V.  = 


M 


P0fhLr 


(4.18) 


where 

t.,  wind  stress 

L relative  characteristic  length  scale 

r for  wind  stress  [=  L/R]  and  L is  a 
length  scale  for  wind  stress  and 
R radius  of  the  earth 
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Substitution  of  these  non-dimensional  variables  into 


(4.10),  (4.11),  (4.4)  and  (4.5)  with  some  rearrangement, 

gives  these  equations  in  terms  of  non-dimensional  variables 
of  the  form 
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where 


Ro  Rossby  number  which  shows  the  relative  importance  of 
local  time  change  term  in  the  equation  of  motion  with 
respect  to  that  of  the  coriolis  term  V. 

1*  *E] 

E horizontal  Ekman  number  which  shows  relative  impor— 

H tance  of  lateral  diffusion  and  coriolis  terms 

Am 

t-  -S,l 

fL 


Ev  vertical  Ekr'in  number,  gives  the  ratio  of  vertical 
diffusion  o.  momentum  to  the  coriolis  term 

.Ai 

fhz 


; . *1 


Pe 


Peclet  number,  shows  importance  of  advection  compared 
to  diffusion  of  density 


V1L 

1= 


Vr  relative  velocity,  ratio  of  barotropic  velocity  to 
baroclinic  velocity 


H relative  depth,  ratio  of  basin  depth  to  thermocline 
r depth 


The  model  is  governed  by  these  six  parameters  assuming  a 
value  for  Lr  is  defined.  All  the  parameters  of  the  resultant 
run  and  constants  of  the  model  are  given  in  table  2. 

Equations  (4.19)  — (4.23)  are  solved  for  variables  u,  v, 
P',  w and  a.  The  variables  u,  v,  and  a are  predicted  from 
(4.19),  (4.20)  and  (4.23).  P'  and  w are  calculated  diagnosti- 
cally from  equations  (4.21)  and  (4.22),  respectively. 
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TABLE  2 . PARAMETERS 

—3 

AND  CONSTANTS 

RO 

4.255  x 10 

Rossby  number 

E„ 

0.532  x 10“2 

Horizontal  Ekman  number 

H 

Ev 

0.851  x 10~2 

Vertical  Ekman  number 

Pe 

10.0 

Peclet  number 

Hr 

50. .0 

Relative  depth 

vr 

0.1 

Relative  velocity 

Ax 

10  x 105  cm 

Zonal  grid  spacing 

Ay 

5 x 105  cm 

Meridional  grid  spacing 

At. 

8 minute 

Time  step 

t 

28.9  day 

Time  scale 

♦ 

O 

O 

Reference  latitude 

n 

2 ir  day 

Rotation  rate  of  the  earth 

R 

6.37  x 108  cm 

Radius  of  the  earth 

g 

—2 

1000  cm  sec 

Gravity  acceleration 

po 

—3 

1.012  gm  cm 

Reference  density 

H 

1 x 105  cm 

Depth  of  the  sea 

L 

10 

Curl  factor  or  relative  length 

r 

scale  for  wind  stress 

K 

. , 2 -1 
1.6  cm  sec 

Eddy  diffusivity 

Kv 

2 —1 
3.2  cm  sec 

Eddy  viscosity 

V1 

4 cm  sec  1 

Characteristic  baroclinic 
velocity 

L 

1 x 107  cm 

Characteristic  length  scale 

am 

0.5  x 108  cm2  sec-1 

Lateral  eddy  viscosity 

7 2—1 

0.4  x 10  cm  sec 

Lateral  eddy  diffusion  coefficient 
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V.  DESCRIPTION  OF  THE  NUMERICAL  MODEL 

A.  SPACE  AND  TIME  DIFFERENCING  TECHNIQUES 

The  space  and  time  differencing  schemes  are  similar  to 
the  schemes  used  by  Haney  (1974)  for  a general  circulation 
model. 

In  the  integration  of  the  primitive  equations  a staggered 
grid  and  centered  space  differencing  scheme  are  utilized. 

The  main  concern  in  choosing  the  staggered  grid  was  avoidance 
of  a computational  mode  which  appears  when  a centered  space 
differencing  scheme  is  used  with  an  unstaggered  grid.  This 
computational  mode  in  space  is  not  present  for  a staggered 
grid.  The  staggered  grid  also  has  advantage  in  saving  comput- 
ing time  due  to  the  presence  of  variables  at  alternate  grid 
points. 

A leapfrog  scheme  is  used  for  time  integration.  The 
starting  scheme  is  a combination  of  forward  and  leapfrog 
scheme  as  shown  below 


Leapfrog  3jAt 


n-n  n=i5 

Forward  ^sAt 


n=l 


(■j  time  step  forward)  + (i  time  step  leapfrog) 


where 


n represents  a time  step 


Time  step  1. 
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In  most  circulation  models  Matsuno  scheme  (Matsuno,  1966) 
is  used  as  a starting  scheme  and  to  remove  the  solution 
separation  which  is  present  when  a leapfrog  scheme  utilized 
exclusively  for  time  integration.  Since  time  integration  in 
this  model  is  not  long  enough  to  cause  solution  separation 
a cheme  such  as  Matsuno  is  not  used  here. 

Each  term  in  the  momentum  and  sigma— t equations  is  evalua- 
ted at  a different  time  step  to  ensure  linear  computational 
stability. 

Neglecting  coefficients,  the  time  differencing  of  the 
equation  of  motion  corresponding  to  (4.20)  is  given  by 


vn+1) 


= v(n  1'  + 2At[PT(n)  + FT(n  1)) 
- 2At  f[Cl  u(n+1)  + C2  u*11-1*] 


(5.1) 


where 

PT  pressure  term 

FT  friction  term,  which  is  composed  of  lateral 

diffusion  of  momentum  and  vertical  eddy  stress 

A trapezoidal  implicit  scheme  is  utilized  for  the  coriolis 

term.  Then,  to  satisfy  consistency  and  linear  computational 

stability,  coefficients  cl  and  C2  must  satisfy  the  relation. 


Cl  + C2  = 1 

j < Cl  < 1 (5.2) 


The  values  Cl  = 0.55  and  C2  = 0.45  are  commonly  accepted 
for  these  constants  which  produce  very  slight  damping  of 
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inertial  oscillations.  Equation  (5.1)  and  an  analogous  equa- 
tion for  u^n+1^  are  solved  simultaneously  for  the  two  unknowns 
u(n+l),  v(n+l)^  T^e  tjJne  differencing  form  Qf  the  sigma— t 
equation  has  the  form 

O (n+l)  = + 2At[ADV^n)+  HD  ^ n”"^ ^ + VD*11-1*] 


♦ «0(cn+l) 


(5.3) 


where 


ADV  advection  term 
HD  horizontal  diffusion  term 
VD  vertical  diffusion  term 


The  leapfrog  scheme  is  conditionally  stable  and  the  time 
step  must  satisfy  the  relation  (for  two  dimensional  wave 
propagation) 


c < 
Ax  — 


(5.4) 


where 


maximum  phase  speed  of  the  waves  present 
in  the  system. 


I it' 


External  gravity  waves  are  removed  and  inertial  oscilla- 
tions are  rendered  neutral  by  (5.1)  and  (5.2).  Therefore 
only  internal  gravity  waves  are  present  to  determine  C.  The 
phase  speed  of  the  internal  gravity  waves  is 


C = /g^H 


(5.5) 


where 

g'  is  modified  acceleration  of  gravity  and 
given  by  the  relation 


where 

g acceleration  of  gravity 

Ap  a typical  value  for  the  vertical  density 

differences 

pQ  reference  density 

—2  —3  —3 

Assuming  g = 1000  cm  sec  , Ap  = 10  x 10  gm  cm  and 

—3  ~ —2 

pQ  = 1.012  gm  cm  gives  a value  for  g'  = 9.88  cm  sec 

Integrations  are  carried  out  with  an  8-minute  time  step  for 

a grid  spacing  of  5 km,  without  any  stability  problem. 

B.  THE  FINITE  DIFFERENCE  EQUATIONS 

The  numerical  method  is  set  down  in  terms  of  the  non- 
dimensional  variables.  The  finite  difference  scheme  is  based 
on  a three-dimensional  array  of  points  with  indices  i,j,k. 
Time  step  is  indicated  by  (n)  as  a superscript.  In  the  nota- 
tion the  letters  i,j,k  are  always  integers. 

Horizontal  spacing  is  such  that  Ax  = 10  km  and  Ay  = 5 km, 
and  uniform  in  the  x and  y directions.  The  horizontal  grid 
pattern  is  shown  in  figure  9.  Where  horizontal  velocity  com- 
ponents (u,v)  are  defined  at  circled  points,  vertical  veloc- 
ity (w)  and  sigma— t (a)  are  defined  at  dot  points. 
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Figure  9.  Placement  of  variables  on 
horizontal  grid  plane.  The  open  cir- 
cles denote  the  definition  points  of 
the  horizontal  velocity  components 
(u,v)  and  dots  denote  definition  points 
of  (a, w) . The  distance  between  adja- 
cent grid  point  is  Ax  = 10  km  and 
Ay  = 5 km  in  x ard  y direction,  res- 
pectively. IM  and  JM  have  values  21 
and  11  respectively. 


For  an  arbitrary  interior  point,  where  (u,v)  are  stored 


longitude  (x)  and  latitude  (y)  are  defined  by 
1 Ax 

. A ri  . i \ i t ^ 


xA  + J ^ U + {i — 1)  ] 


j + 7 - ¥ 11  + 


i = 1,2,...,  IM— 1 


j = 1,2,...,  JM— 1 


(5.7) 


The  vertical  pattern  of  variables  is  shown  in  figure  10. 
The  vertical  index  indicated  with  k and  the  vertical  velocity 
w are  given  along  the  surface  and  bottom  and  along  the  layer 
boundaries.  The  variables  (u,v,a)  are  located  at  the  kth 
level.  Depth  of  each  level  is  defined  as 
AZ. 


Z1  =--t 


Zk  - Z1  - £ (AZm  + AZjn_1 ) /2  k = 2, 3...  KM 


(5.8) 


Layer  thicknesses  for  each  level  are  given  in  table  3. 
Table  3.  LAYER  THICKNESSES 


k 

1 2 3 4 5 

6 

7 

8 

9 

AZk 

— 

0.25  0.25  0.3125  0.6875  1.0 

1.5 

6.0 

15.0 

25.0 

To  resolve  the  surface  layer  more  accurately,  four  levels 
are  located  within  20  meters.  Close  to  the  bottom  a coarser 
grid  is  used.  Resolution  in  the  surface  layer  is  greater  than 
in  the  bottom  layer. 


#» 


6 U,V,T 6 0 0 

6 Vj _ w_  _ 

Ksa?  u-y.T !0  0 0 


< 1 Figure  10.  Vertical  structure  of  the 

model  (u,v,a)  are  defined  at  integer  K 
values  and  w is  defined  along  the  layer 
f 1 boundaries. 

H 


;/  " < 


The  finite  difference  form  of  the  hydrostatic  equation 
has  the  form 


P . . , = (a— a ) . . , AZ, 

1]  1 o h 


Pijk  “ Pij  , k— 1 + (a_co)  ijk-^s  *AZk-is  k_2'3' 


KM 


(5.9) 


and  the  vertical  average  of  the  departure  pressure  is 
calculated  with  the  relation 

k=KM 

P.  . = V P . AZ, 

n *-•  13k  k 


(5.10) 


k=l 


The  departure  P'  is  calculated  from 


P!  = P.  - P.  . 
ljk  ljk  lj 


(5.11) 


which  is  the  finite  difference  analog  of  (4.23). 

The  finite  difference  approximation  of  the  continuity 
equation,  in  which  the  vertical  velocity  between  the  levels 
is  calculated,  has  the  form 


w.  = w.  ..  1 + [U  + V ] . AZ, 

ijk+h  13k— x y i]k  k 


where 


and 


Ux  ijk  55  ^ux  j+h  + Ux  j— ^ ik 


Vy  ijk  = ls(vy  i+H  + ^ 


ux  ij+^k  Ax^ui+%  ui— ^j+3gk 


ry  i+^sjk  Iy(vj+Js  vj-^)i+isk 


(5.12) 


(5.13) 


(5.14) 
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which  are  analogs  of  and  (|^)  respectively.  Gradients 

of  horizontal  velocity  components  u and  v are  defined  at  half— 

x y 

integer  grid  points  west  and  south  of  the  u and  v storage 
points  shown  in  figure  11a. 

Integration  from  the  surface  downward,  using  surface 
boundary  condition  w=0,  is  done  by 


wi jk+% 


- r <ux + Ytj Az 


(5.15) 


The  formula  for  the  finite  difference  approximation  of 
the  equation  (4.19)  by  which  u is  predicted  is 


+ 


+ 


n+1  n— 1 

ru  ~ u i 

1 2At  1 i+hj+hk 

P • p*  P 

i+1  j Pijin  + 

Ax  1 k Ro 


uy  j+h  j + 1 Uy  i+^sj 

Ay 


p ' — p ' 

1 r 1+1  j+1  1 ]+l 

2Ro  1 Ax 


rUX  i+1  j+^s  ux  i 

L Ay 

, n— 1 fv  ruz  k—h  ~ uz  k+h  n-1 

k Ro  AZ^  i+h  j+h 


+ ife  tCl  v°+1  + C2  \+%  j+%k 


(5.16) 


where 

and 


uyj+l 

uzk— % 


ui+ih  ui+h) 

Ay  ' i+^k 

uk— 1 ~ uk' 

AZk_^  ' i+hj+h 


(5.17) 

(5.18) 


Equations  (5.16)  and  (5.17)  represent  (|^-)  and  (|j)  and 
are  located  half  grid  distance  south  and  above  the  velocity 
(u)  storage  points,  respectively,  as  shown  in  figure  11  a— b. 


The  formula  for  the  finite  difference  approximation  of 
the  equation  (4.20)  in  which  v is  predicted  has  the  form 


n+1  n — 1 p < _ pi 

rX 1 - _ 1 r i+1  i+1  P i+li 

2At  Ji+Jsj+35k  2Ro  1 

+ " Pijjn  + fg  cvx  i+1  j+%  ~ vx  ii+h 


Ay 


Ro 


Ax 


+ VY  i+^j+l  vy  i+%jT  (n— l) 
Ay  J k 


E 

+ ■ 


Vz  k— vz  k+% 
AZ, 


i+H>j+% 


+ Ro  ST  tY  ” ik  [Clun+1  + C2un  X]  . 

Ro  Hr  i+Jj  Ro  1 i+hj+hk 


where 


Vx  i+1 
Vz  k— % 


vi+lh  vi+k] 

Ax  j+^k 

Vk— 1 ~ vk. 

Azk -3s  i+?23+}5 


(5.19) 


(5.20) 

(5.21) 


Equations  (5.20)  and  (5.21)  are  analogues  of  (|X)  and 
and  are  stored  south  and  above  the  velocity  storage 
points,  as  shown  in  figure  lla-b. 
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The  finite  difference  approximation  of  the  sigma— t equa 
tion  has  the  following  form. 


n+1  n— 1 

a — o 


(n)  . 1 ,°x  i+^j+ig  ax  i— %j+h 


+ y i+ki+h  y i+^j-^n-l 
Ay  7 k 


_K 1_  _ ,q2  k+^s  az  k— n— 1 

K Ro  vl  AZ,  ' i+hj+h 

V Jv 


+ 5c(qn+1)ijk 


Aijk  4Ax  ^ui+Jgj+is  + Ui+}sj— ^ h1  ^aij  + ai+lj^k 


^ui— hj+h  + ui— igj— ’ ^ai— lj  + aij^k^ 


+ -t-t — [ (v.,,.,,  + v._i  . , ) • (o  . . , + a . . ) . 

4Ay  l+^j+^s  i—hj—h  lj— 1 ij  k 


— (v.  , , . , + v.  ,)  • (a  • . , +a..),] 
i+hj—%  i— V i]— 1 i]  k 


+ 7AZ.  ^ (wijk-*^  * (aijk  + aijk— 1^ 


^wijk+^  ^°ijk  + °ijk+l^  ^ 


ax  i+%jk  - Ax  ^ai+l  ai^jk 

°Y  ij+^k  - 'Ey  ^aj+l  aj^ik  (5.24) 

az  ijk+^s  ~ TTzJ ^azk  azk+l^  ij 

As  shown  in  figure  (11  a— b)  ax  and  and 

are  defined  at  half  grid  distances  south,  east  and  above  the 
sigma— t storage  points. 

Before  proceeding  to  a new  time  step  convective  adjustment 
is  applied.  It  is  difficult  to  state  the  convective  adjust- 
ment mechanism.  With  this  mechanism  the  density  structure  is 
forced  to  remain  stable.  It  may  be  expressed  by 

0 °Z  < 0 

« = Q Z (5.25) 

°Z  > ° 

In  case  of  unstable  stratification,  vertical  mixing  be- 
comes effectively  infinite.  At  each  -time  step  this  infinite 
mixing  is  included  in  numerical  solution  by  testing  sigma— t 
profile  for  unstable  lapse  rate  (o^^  > cr^)  . If  this  condi- 
tion exists  new  values  of  sigma— t for  those  two  layers  are 
set  equal  to  the  vertical  average  sigma— t of  the  two  layers 
instantaneously.  This  process  is  repeated  until  complete  sta- 
bility is  reached  for  the  entire  layer. 


C.  SPECIFICATION  OF  BO UN DAT I ON  CONDITIONS 


As  shown  in  Figure  3—1  o,w  are  variables  which  are  de- 
fined at  lateral  boundaries;  also  horizontal  gradients  of 
(u,v)  are  defined  at  these  lateral  boundaries.  The  boundary 
conditions  are  satisfied  by  the  following  relations  which 
are  given  for  eastern  and  western  boundaries.  Similar  con- 
ditions exist  for  north  and  south  boundaries  with  the  excep- 
tion of  open  boundaries. 

Normal  components  are  set  equal  to  zero  by 


u 


i=l 


= iH.) 

Ax' i=% 


u 


X 


i=IM 


= •^) 

Ax' i=IM-% 


(5.26) 


In  the  momentum  equation  zero— slip  condition  is  applied 
to  the  valocity  tangent  to  the  boundary.  On  the  other  hand, 
when  computing  the  advection  term  in  the  density  equation 
free  slip  condition  exists  for  the  velocity  tangent  to  the 
boundary.  This  implies  the  following  conditions 


u 


i=l 


u 


i=lh 


(5.27) 


Ui=IM  ~ UIM— % 
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Definition  of  vertical  gradients  of  (u,v)  allows  the 
writing  of  the  vertical  stress  term  in  flux  form  easily. 
Boundary  conditions  at  the  surface  and  bottom  are  defined 


in  the  following  manner: 


uz  k - f - 0 


(5.28) 


Zonal  wind  stress  is  assumed  zero  and  only  a constant 


meridional  wind  stress  is  defined. 

. V L 

v k = ± = ,r  r T y 

v r?  ^ ^ L _ 


(5.29) 


On  the  other  hand  it  is  assumed  that  horizontal  velocities 
vanish  at  the  bottom.  The  following  conditions  force  veloc- 
ities (u,v)  to  reach  zero  at  the  bottom 


Z k=KM+% 


2u  | 

AZ 1 k=KM 


(5.30) 


Z k=KM+?s 


— I 

AZ 1 k=KM 


The  boundary  conditions  of  zero  mass  fluxes  across  the 
side  walls  and  bottom  are  satisfied  by  imposing  the  following 
conditions : 


(0u)i=!i  " -«’u>i.3/21jk 


(0“)  i=iH+lj  ' (au)  jk 


(5.31) 


Similar  conditions  exist  at  north  and  south  boundaries. 

In  the  diffusion  term  zero  flux  and  insulation  are 
satisfied  by  the  following  conditions: 


— a 


i*3/2 


(5.32) 


o = — a 

Xi=IM+^  Xi=IM— ^ 

where  ax  is  defined  at  half  grid  distances  west  of  the 
sigma— t storage  points,  as  shown  in  figure  11a.  Similar 
conditions  exist  at  north  and  south  for  meridional  gradients 
of  sigma— t,  with  the  exception  of  open  boundaries  where 
sigma— t is  defined  outside  the  domain. 

At  the  sea  surface  downward  flux  is  absent. 


az  k = 


(5.33) 


The  sigma— t equation  is  solved  for  only  the  upper  six 

levels.  It  is  assumed  that  sigma— t for  the  rest  of  the  depth 

range  has  a value  o , which  is  held  constant  during  integra— 

c 

tion.  This  enters  the  system  in  the  calculation  of  the 
vertical  gradient  of  sigma— t one  half  grid  distance  below 
the  sixth  level. 


a 


Z 


k=6*5 


k=6 


and  ac  has  a value  28.5  in  sigma— t units. 


(5.34) 
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VI . RESULTS 


Computations  are  carried  out  for  different  values  of  the 
parameters  with  an  integration  period  of  24  hours.  Parameters 
of  these  runs  are  given  in  Table  2.  Results  of  the  model  are 
given  for  the  fields  of  horizontal  velocity  components  (u,v) , 
vertical  velocity  (w)  and  density  (sigma— t)  in  dimensional 
form. 

The  horizontal  velocity  field  at  z = — 2.5  m is  given  in 
Figure  12.  The  horizontal  stress  exerted  on  the  sea  surface 
by  southerly  wind  causes  the  surface  water  to  move  in  a 
generally  north— east  direction.  This  is  a consequence  of 
movement  of  surface  water  as  an  Ekman  layer,  drifting  to  the 
right  of  the  wind  stress  in  the  Northern  Hemisphere.  Since 
the  curl  of  the  wind  is  zero  divergence  of  the  Ekman  drift  is 
also  zero.  But  horizontal  velocity  is  strongly  convergent  and 
divergent  near  the  boundaries.  Northeastward  drift  takes 
place  over  the  entire  basin  at  this  depth.  The  strongest  flow 
has  a magnitude  20  cm/sec  in  the  central  part  of  the  basin  and 
diminishes  toward  the  boundaries.  Near  open  boundaries 
velocities  are  smaller  compared  to  the  velocities  at  neighbor^ 
ing  grid  points.  The  southward  flow  at  these  points  repre- 
sents the  exchange  with  the  Black  Sea  water  at  the  north  and 
Mediterranean  water  at  the  south.  The  magnitudes  of  these 
currents  are  5 cm/sec  and  2 cm/sec  at  grid  points  just  near 
the  entrances  of  the  Straits  of  Canakkale  and  the  Strait  of 
Istanbul.  The  applied  meridional  wind  stress  has  a value  of 
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0.3  dyn  cm  which  is  not  strong  enough  to  move  surface  water 
to  the  north  as  it  does  at  neighboring  grid  points. 

In  the  second  layer  z = —7.5  m the  horizontal  velocity 
changes  direction  to  the  right  and  becomes  weaker.  Horizontal 
velocities  at  that  level  are  shown  in  figure  13w;  they  have  a 
magnitude  5 cm  sec  1.  Flow  adjacent  to  the  straits  is  still 
southerly  and  has  a magnitude  5 cm  sec  ^ and  3 cm  sec  ^ at 
the  vicinities  of  the  Straits  of  Istanbul  and  Canakkale 
respectively. 

The  circulation  pattern  at  lower  levels  differs  from  the 
upper  level  flows  both  in  magnitude  and  direction.  These 
differences  can  be  observed  in  the  circulation  pattern  at 
the  z = —40.0  m level  (Figure  14).  The  pattern  is  not  irregu- 
lar compared  to  the  upper  layers. 

The  northeastward  surface  drift  in  the  first  layers  is 
accompanied  by  a westward  and  partly  southwestward  flow  in 
the  bottom  layers.  A northward  flow  adjacent  to  the  Strait 
of  Canakkale  brings  Mediterranean  water  into  the  basin.  A 
similar  northward  flow  of  Mediterranean  water  does  not  appear 
at  the  northern  strait  due  to  the  shallow  sill  depth  at  that 
point. 

Flow  in  the  layers  close  to  the  bottom  becomes  weaker. 

The  influence  of  the  circulation  at  the  straits  on  the  general 
pattern  is  small.  The  circulation  pattern  at  a depth  700.0  m 
is  shown  in  figure  15.  The  maximum  magnitude  is  0.033  cm  sec 

One  advantage  of  the  numerical  modeling  study  is  that 
it  gives  estimates  of  important  oceanographic  variables  that 
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cannot  be  obtained  by  measurements.  One  such  variable  is 
the  vertical  velocity.  The  vertical  velocity  at  the  base 
of  the  first  layer  is  given  in  Figure  16.  The  magnitude  of 
the  vertical  velocity  differs  markedly  between  the  boundaries 
and  the  interior.  Rising  motion  (upwelling)  is  taking  place 
in  the  south  and  southwest  and  sinking  motion  (downwelling) 
is  taking  place  in  the  north— northeast.  Vertical  velocity 
is  mainly  determined  by  the  Ekman  drift  current  intersecting 
the  boundaries.  Southerly  wind  stress  produces  upwelling 
and  downwelling  at  the  south  and  north,  respectively.  Up- 
welling  extends  to  the  north  on  the  western  side  and  down- 
welling  to  the  south  on  the  eastern  side.  The  maximum  magni- 
tude of  the  vertical  motion  is  358  cm  day  1 . Strong  upwell— 
ing  and  downwelling  take  place  adjacent  to  the  straits.  This 
is  consistent  with  the  nature  of  the  horizontal  current 
system  in  the  vicinity  of  the  straits. 

Vertical  velocity  at  the  bottom  of  the  second  layer, 
z = — 7.5  m,  is  shown  in  Figure  17.  The  general  pattern  is 
similar  to  the  vertical  velocity  at  the  base  of  the  first 
layer.  Strong  vertical  velocity  is  present  along  the  lateral 
walls  to  the  southwest  and  northeast  where  it  has  a magnitude 
+ 300  cm  day""1.  With  depth  the  effect  of  wind  on  vertical 
velocity  decreases  due  to  the  vertical  stability.  The  verti- 
cal velocity  pattern  at  the  base  of  the  eighth  level  is 
shown  in  Figure  18. 

The  sigma— t pattern  at  the  first  level,  2.5  m is  shown 
in  Figure  19.  This  pattern  at  z = — 2.5  m is  specified  by 
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Figure  17.  Vertical  velocity  at  the  base  of  the  second 
layer  (10.0  meter). 
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the  surface  boundary  condition.  At  this  level  the  surface 
sigma— t is  a function  of  latitude  only  and  held  constant 
during  the  integration  period. 

Deviations  from  the  surface  pattern  begin  to  appear  just 
below  the  first  level  at  7.5m.  The  sigma— t pattern  for  the 
second  level  is  shown  in  Figure  20.  A close  relationship 
is  present  between  this  sigma— t pattern  and  the  vertical 
velocity  field  shown  in  Figure  17.  At  the  north  and  northeast 
sinking  motion  brings  in  low  density  and  water  and  at  the 
south  rising  motion  brings  in  high  density  water  to  this 
depth.  This  pattern  is  more  strongly  tied  to  the  surface 
boundary  conditions  on  sigma— t than  at  other  levels.  There 
is  a low  density  water  pool  at  the  north.  As  a consequence 
of  upwelling  a high  dense  water  pool  is  present  at  the  south 
extending  in  an  east— west  direction.  Magnitudes  of  sigma— t 
are  16.5  and  18.0  at  north  and  south  respectively. 

The  sigma— t pattern  at  the  fourth  level  z = — 20  m is 
shown  in  Figure  21.  The  effect  of  vertical  motion  appears 
as  a low  density  water  pool  at  the  northeast.  Below  this 
level  z = —40  (not  shown),  sigma— t is  uniform  and  has  a 
value  28.5.  This  is  consistent  with  the  sigma— t pattern 
given  in  Figure  6. 

Further  insight  into  these  results  may  be  gained  by  examin- 
ing the  meridional  and  east-west  cross  sections  of  sigma— t 
and  the  (u,v)  components  of  the  horizontal  velocity.  There 
is  also  the  possibility  of  comparing  the  predicted  sigma— t 
pattern  to  the  observations. 
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Figure  20.  Horizontal  sigma— t pattern  at  second  level 
(7.5  meter) . 


(20  meter) 


Meridional  cross  sections  of  sigma— t at  different  longi- 
tudes are  shown  in  Figures  22—24.  The  general  patterns  are 
similar.  High  density  water  rises  to  the  surface  and  strati- 
fication becomes  shallower  further  south.  Below  40  m density 
is  constant.  This  is  consistent  with  the  observed  sigma— t 
field.  At  the  west  and  south  very  strong  upwelling  and  weak 
downwelling  at  north  are  shown  in  Figure  22.  In  the  central 
region  there  is  upwelling  in  the  south  and  downwelling  in  the 
north  as  shown  in  Figure  23.  Further  east  the  effect  of 
strong  downwelling  is  shown  by  the  rate  of  tilting  of  the 
isopycnals  at  the  north  compared  to  the  south  in  Figure  24. 

Meridional  cross  sections  of  zonal  velocity  in  the  east, 
central  and  western  parts  of  the  basin  are  shown  in  Figures 
25—27.  The  zonal  flow  at  the  surface  is  eastward  and  changes 
direction  with  depth.  Zonal  flow  is  westward  at  greater 
depths.  Strong  flow  is  present  at  the  surface  and  decre?ses 
rapidly  with  depth.  Variation  with  longitude  can  be  observed 
by  comparing  the  figures.  In  the  central  part  zonal  velocity 
has  a magnitude  16.0  cm/sec  at  the  surface  and  decreases 
both  laterally  and  vertically.  Below  20  meters  flow  is  west- 
ward and  has  a magnitude  0.6  cm  sec  \ 

East— west  cross  sections  of  sigma— t patterns  for  three 
latitudes  are  shown  in  Figures  28—30.  A regular  stratifica- 
tion is  deformed  by  upwelling  and  downwelling  of  the  eastern 
and  western  sides  of  the  basin  (Figure  29) . Departures  are 
present  in  the  north  and  south,  as  shown  in  Figures  28  and  29 
respectively.  Intrusion  of  less  dense  Black  Sea  water  can 
be  identified  in  Figure  30. 
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Figure  22.  Meridional  cross  section  for  sigma— t at  the 
western  part  of  the  basin. 


The  vertical  cross  sections  of  the  meridional  velocity 
which  correspond  to  the  previous  sigma— t sections,  are  given 
in  Figures  31—33.  There  is  a northward  flow  at  the  surface 
in  each  section.  The  cross  section  near  the  southern  bound- 
ary shows  a southerly  flow  below  the  surface  layer.  Flow 
changes  direction  at  about  20  meters.  In  the  vicintiy  of 
the  Strait  of  Istanbul  the  effect  of  the  dpen  boundary  is 
indicated  by  a weak  northerly  flow.  There  is  a southerly 
flow  below  that  and  at  depth  about  30  m another  reversal 
takes  place.  Weak  northward  flow  is  associated  with  the 
southerly  wind.  The  southward  and  northward  flows  at  depth 
are  consistent  with  observations. 


Figure  32.  East— west  cross  section  of  the  meridional 
velocity,  (v) , at  the  central  part  of  the 


VII.  CONCLUSIONS 


In  order  to  investigate  the  dynamics  of  the  Sea  of  Mar- 
mara a .circulation  model  is  formulated  based  on  the  hydro- 
static and  Boussinesq  approximations.  Since  very  little  is 
known  about  the  turbulence  characteristics  of  the  ocean  or 
adjacent  seas,  viscosity  and  conductivity  are  replaced  by 
new  terms  representing  the  contributions  of  the  smaller  scale 
motions  to  the  exchange  of  momentum  and  density.  Horizontal 
eddy  transfer  of  momentum  is  accomplished  using  a constant 
eddy  viscosity  and  diffusivity  coefficient.  For  the  exchange 
of  momentum  and  density  in  the  vertical  a constant  vertical 
diffusion  coefficient  also  is  used.  The  vertical  mixing  of 
density  is  enhanced  by  the  use  of  an  instantaneous  adjustment 
mechanism  to  eliminate  unstable  lapse  rates.  A lateral  mix- 
ing coefficient  for  density  is  less  than  that  used  for 
momentum.  These  coefficients  should  be  sufficiently  small 
so  as  not  to  obscure  lateral  transfer  of  the  quantities,  and 
a non  zero  momentum  coefficient  is  needed  to  satisfy  the 
zero  slip  boundary  condition.  In  this  study  constant  values 
of  the  coefficients  are 

7 2—1 

Ay  = .4  x 10  cm  sec  (density) 

A 7 2—1 

M = 5 x 10  cm  sec  (momentum) 


and 


The  vertical  diffusion  coefficient  which  is  difficult  to 

evaluate  either  by  measurement  or  in  principle,  is  assigned 
2 —1 

1.6  and  3.2  cm  sec  for  density  and  momentum  respectively. 
Other  parameters  and  constants  are  given  in  Table  II. 

Further  study  of  thermocline  penetration,  water  exchange, 
vertical  mixing  processes  and  current  gradients  might  provide 
a better  possibility  for  matching  these  coefficients. 

Although  many  important  factors,  such  as  the  irregular 
bottom  topography,  bottom  friction  and  non— linear  terms  in 
the  equation  of  motion  have  been  omitted,  the  results  from 
the  model  are  generally  encouraging.  More  realistic  results 
may  be  achieved  by  incorporating  these  effects  systematically 
in  further  studies. 

>'s 

— - p Integration  in  the  model  is  carried  out  with  a constant 
southerly  wind  stress  which  was  chosen  arbitrarily.  Although 
the  scale  of  the  motion  is  small  compared  to  a characteristic 
length  scale  for  meteorological  features,  wind  stress  may  be 
significantly  variable  over  the  sea  due  to  the  geographic 
location  of  the  sea.  Also  seasonal  variations  of  the  wind 
may  have  a considerable  effect  on  the  vertical  mean  part  of 
the  current  which  is  neglected  during  this  study. 

The  model  does  not  separately  predict  temperature  and 
salinity  which  are  sometimes  important  quantities.  Diagnostic 
calculation  of  the  density  with  the  aid  of  predicted  tempera- 
ture and  salinity  fields  can  be  done  with  a very  simple 
modification  of  the  model.  This  would  require  a calculation 
of  the  heat  and  salt  fluxes  at  the  surface  and  meridional 


1 


salt  and  heat  separately  in  the  model.  It  would  be  more 
difficult,  however,  to  define  boundary  conditions  for  these 
quantities  at  the  open  boundaries.  It  is  also  difficult  to 
simulate  seasonal  variations  of  these  quantities  at  the 
open  boundaries.  Once  these  are  defined  properly,  based 
on  accurate  observations,  there  are  no  inherent  difficulties 

I 

in  simulating  the  dynamical  processes  in  the  sea  using  the 

general  principles  involved  in  this  model.  Even  though  this 

model  depends  on  the  particular  hypothesis  used  for  hori— 

i 

zontal  and  vertical  eddy  transport  of  heat,  salt  and  momentum, 
?'  it  is  clear  that  valuable  studies  on  small  scale  water  bodies 


APPENDIX  A 


COMPUTER  PROGRAM  DESCRIPTION 

The  computer  program  is  written  in  FORTRAN  IV  and  was 
used  with  the  IBM  360/67  computer  system  at  the  W.  R.  Church 
Computer  Center,  Naval  Postgraduate  School. 

The  overall  program  is  divided  into  two  basic  subpro- 
grams (1)  the  main  program  and  associated  subroutines  (2)  an 
access  program  which  draws  and  writes  the  results  of  the 
first  program. 

The  main  program  consists  of  nine  subroutines  which 
calculate  different  terms  of  the  equation  of  motion  and  the 
density  equation,  pressure,  vertical  velocity  and  changes 
variables  for  next  time  step. 

FORTRAN  IV  symbols  for  the  primary  program  and  a brief 
description  of  the  subroutines  are  given  below: 

UMI  Zonal  velocity  component  at  n— 1 time  step 

U Zonal  velocity  component  at  n time  step 

UA1  Zonal  velocity  component  at  n+1  time  step 


VM1 

Meridional 

velocity 

component 

at 

n— 1 

time 

step 

U 

Meridional 

velocity 

component 

at 

n 

time 

step 

UAl 

Meridional 

velocity 

component 

at 

n+1 

time 

step 

SGM1 

Sigma— t at 

n— 1 time 

step 

SGMT 

Sigma— t at 

n time 

step 

SGMT 

Sigma— t at 

n+1  time 

step 

ADV 

ADV 

ADV 

A 

PBAR 

W 

UXG 

UYG 

UZG 

VXG 

VYG 

VZG 

SGX 

SGY 

SGZ 

TX 

TY 

DZ 

DZ 

BB 

CA 

BO 

CD 

DX 

DY 

RO 


Total  advection  term  in  density  equation 

Local  rate  of  change  of  sigma— t 

Pressure 

Surface  pressure 

Vertical  average  pressure 

Vertical  velocity 

Gradient  of  zonal  velocity  in  x direction 

Gradient  of  zonal  velocity  in  y direction 

Gradient  of  zonal  velocity  in  vertical 

Gradient  of  meridional  velocity  in  x direction 

Gradient  of  meridional  velocity  in  y direction 

Gradient  of  meridional  velocity  in  vertical 

Gradient  of  sigma— t in  x direction 

Gradient  of  sigma— t in  y direction 

Gradient  of  sigma— t in  vertical 

Meridional  wind  stress 

Zonal  wind  stress 

Layer  depth 

Layer 

Depth  dependent  velocity  at  the  Northern 
fictitious  boundary 

Depth  dependent  velocity  at  the  Southern 
fictitious  boundary 
fictitious  boundary 

Sigma— t defined  outside  the  domain  at  North 
Sigma— t defined  outside  the  domain  at  South 
Horizontal  grid  spacing  in  x direction 
Horizontal  grid  spacing  in  y direction 
Rossby  number 


EH 

EV 

PEI 

W1 

AK 

AKV 

RL 

HH 

DT 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 


Horizontal  Ekman  number 
Vertical  Ekman  number 
P6clet  number 
Relative  velocity 
Vertical  eddy  diffusivity 
Vertical  eddy  viscosity 
Curl  factor 
Depth  of  the  basin 
Time  step 

PRES  Calculates  pressure 

VERW  Calculates  vertical  velocity 

ADVEC  Calculates  advection  term  in  the 

density  equation 

HDGRD  Calculates  gradients  of  the  sigma— t 

DIFFO  Calculates  local  time  rate  of  the 

change  sigma— t by  subtracting  advec— 
tion  term  from  the  diffusion  term 

SIGEQ  Calculates  new  sigma— t and  makes 

convective  adjustment  for  unstable 
lapse  rates. 

HVGRD  Calculates  gradients  of  the  hori- 

zontal velocities 

HORV  Calculates  horizontal  velocity 

components,  u,v. 


A descriptive  flow  diagram  of  the  program  is  shown  in 
Figure  34  . 
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